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Abstract. 

A lattice model of three-state stochastic phase-coupled oscillators has been shown 
by Wood et al. (2006 Phys. Rev. Lett. 96 145701) to exhibit a phase transition at 
a critical value of the coupling parameter, leading to stable global oscillations. We 
show that, in the complete graph version of the model, upon further increase in the 
coupling, the average frequency of collective oscillations decreases until an infinite- 
period (IP) phase transition occurs, at which point collective oscillations cease. Above 
this second critical point, a macroscopic fraction of the oscillators spend most of 
the time in one of the three states, yielding a prototypical nonequilibrium example 
(without an equilibrium counterpart) in which discrete rotational (C3) symmetry is 
spontaneously broken, in the absence of any absorbing state. Simulation results and 
nucleation arguments strongly suggest that the IP phase transition does not occur on 
finite-dimensional lattices with short-range interactions. 
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1. Introduction 



Models of coupled oscillators exhibit a surprising range of symmetry-breaking transitions 
to a globally synchronized state. In the paradigmatic Kuramoto model, for instance, 
oscillators with intrinsic frequencies Ui coupled via their continuous phases 9i can exhibit 
stable collective oscillations, with breaking of time translation symmetry [U Ej [3l IH [5] . 
When Ui = 0, the Kuramoto model acquires an additional 6i — —6i reflection symmetry, 
which can also be broken at the onset of synchronization, yielding surprising response 
properties [6]. 

Within the class of discrete-phase models, the paper-scissors-stone game is an 
example of a system with three absorbing states which can present either global 
oscillations or spontaneous breaking of discrete rotational (C3) symmetry [71 [8| [9| [TOt [TT] . 
More recently. Wood, Van den Broeck, Kawai and Lindenberg proposed a family of 
models of phase-coupled three-state stochastic oscillators which can undergo a phase 
transition to a synchronized state [121 lEl [HI HH] for sufficiently strong coupling. From 
here on each of them is called Wood et a/.'s cyclic model (WCM). Although the WCM 
also has C3 symmetry, it has no absorbing state. 

As noted previously [H], the first WCM [12] has an unusual behavior after the 
phase transition to a synchronized state. As the phase coupling is further increased, 
collective oscillations slow down while oscillators remain partially synchronized. Here 
we have a closer look at this phenomenon, in particular addressing whether or not there 
is a second phase transition with breaking of C3 symmetry. 

The paper is organized as follows. In section [2] we review the WCM and the 
essentials of the known transition to a synchronized phase. Our results regarding a 
putative second phase transition are described in section |3l We conclude in section jH 
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2. Model 

In the WCM, the phase (px at site x {x = 1, . . . , A^) can take one of three values: 0^ = 
j3;27r/3, where G {0, 1,2} (for simplicity, from here on we employ x to denote sites 
and i to denote states — always modulo 3). The only allowed transitions are those 
from i to j + 1 (see Fig. [1]), with transition rates 

'a (iV,-+i - N,) 



9j,j+i = 9 exp 



(1) 



where g is a. constant (which can be set to unity without loss of generality), a is the 
coupling parameter, Nj is the number of nearest neighbors in state j, and z is the 
number of nearest neighbors. The rates are invariant under cyclic permutation of the 
indices, i.e., belong to the symmetry group C3 of discrete rotations, and strongly violate 
the detailed balance condition. 

The mean-field approximation is obtained by replacing Nj/z in the argument of 
the exponential of Eq. [T] with the corresponding probability, Pj , leading to the following 
set of equations: 

= 9j~i,jPj-i - 9j,j+iPj > (2) 

where now gjj+i = gexp[a{Pj+i — Pj)]. These also represent the equations for a 
complete graph in the limit N 00: since z = N — 1, we can replace Nj/N with Pj in 



1 







g 







Figure 1. Transition rates for an isolated unit. 




Figure 2. (Color online) Phase space {Pi,P2). Grey triangles correspond to the 
allowed region < P2 < 1 ^ ^'i- a increases from left to right. In the upper 
panels, dashed (solid) lines denote the nuUclines Pi = {P2 — 0). Lower panels show 
trajectories from complete graph simulations with N — 5000 (thick lines) and solutions 
of the mean- field equations (thin lines with arrows). Panel (a) shows that Pjyg is the 

only fixed point in the system for a = 1.6 > Oc. (b) Pi/^ is unstable, with trajectories 
converging to a stable limit cycle, which corresponds to sustained global oscillations 
(see also Fig.[3K)- (c) For a = 3, the nullclines approach one another in three additional 
positions, thereby creating three symmetric saddle-node ghosts, (d) The limit cycle 
is almost as large as the perimeter of the allowed region. Note the overcrowding of 
lines near the ghosts, with consequent anharmonicity in the time series (Fig. [3})). 
(e) For a = 3.25 there are six additional fixed points which are born simultaneously 
in three saddle-node bifurcations at a = a^. (f) In this case, if the system starts from 
the symmetric unstable fixed point -Pjyg, fluctuations determine to which stable fixed 
point it goes. Closed (open) symbols denote stable (unstable) fixed points. 

the infinite-size limit. Normalization reduces these equations to a two-dimensional flow 
in the (Pi, P2) plane (Fig. E]). 

Starting from the complex order parameter, 




x=l 



one can define quantities that characterize the collective behavior of the system. A 
configuration with v ^ has unequal numbers of sites in the three states, which led 
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Figure 3. (Color online) Panels (a) and (b) show the evolution of Pi for a = 1.6 
and a = 3, respectively. Points represent simulations for complete graph with N = 
10000 whereas lines show the mean-field solution, (c) Dependence of r and w on a, 
exhibiting the two phase transitions, (d) Xr versus a showing peaks at the transitions 
[system sizes as in (c)]. 
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Kuramoto to propose 

r^m)t)s (4) 
as the order parameter for synchronization [I],[3l[T2], where 0^ denotes a time average 
over a single reahzation (in the stationary state), and 0^ an average over independent 
reahzations. Note that r > is consistent with, but does not necessarily imply, globally 
synchronized oscillations. The latter is characterized by a periodically varying phase 
of [161[I71[I81[T9]. 

In the mean- field equations ([2]), the transition to the synchronized regime is 
associated with a supercritical Hopf bifurcation at a = Oc = 1.5. The trivial fixed 
point P*/3 = (Pi = 1/3, P2 = 1/3) loses stability at a = Oc, and a limit cycle 
encircling this point appears (Fig. [2|i-b) [I2]. For a > Oc, sustained oscillations 
in Pj characterize synchronization among the oscillators (Fig. |3^). Correspondingly, 
r grows continuously ~ (a — a^l^ at the transition (Fig. [31:;), with a mean-field 
exponent (3 = 1/2 [12]. The scaled variance Xr = L'^ K'^I^Dt)* ~ diverges with the 
system size at criticality, as shown for simulations of complete graphs using L*^ = in 
Fig. Hi. 

Note that the phase transition to stable collective oscillations is associated with 
the breaking of the continuous time-translation symmetry: the time series Po(t), Pi(^) 
and P2(t) become periodic for a > a^. But since they are statistically identical, except 
for a phase, C3 symmetry still holds. Increasing a above Oc enhances synchronization 
among the oscillators, leading to increasing oscillation amplitudes, as shown in Figs. [2ji 
and |3b. 

3. The Second Phase Transition 

3.1. Mean- field theory and the complete graph 

Wood et al. found that the increasing amplitude of oscillation is accompanied by 
a decreasing angular frequency u = 2tt/{t), where (r) is the mean time between 
peaks in P^ (Figs. |3K-c). This can be understood qualitatively from the exponential 
dependence of the rates of Eq. [1] on the neighbor fractions: when a state is highly 
populated, the rate at which oscillators leave this state becomes very small. In addition, 
the rate at which a highly populated state attracts oscillators from its predecessor is 
also very high. 

We study the effect of a further increase in the coupling a, beyond the range of 
values investigated in Ref. [I2], using mean- field theory and analysis on a complete 
graph of sites. Note that the state of the latter is completely specified by two integer 
random variables, A^^i and A^2- We study the complete graph (for up to 1000) via exact 
(numerical) stationary solution of the master equation [20] and (for ranging from 10 
to 10^) via event-driven Monte Carlo simulations. The results of the two approaches 
are fully consistent. 
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In mean-field theory, when a reaches an upper critical value a'^ ~ 3.102 439 915 64, 
three symmetric saddle-node bifurcations occur simultaneously, and the period of the 
collective oscillations diverges. Above a^, there are three symmetric attractors in 
the system, and 3-fold rotational (C3) symmetry is spontaneously broken: Figure |2]F 
illustrates how the final fate of three different trajectories starting from the same initial 
conditions is determined by fluctuations. Moreover, the transition can be considered 
reentrant, in the sense that time-translation symmetry, which had been broken at Oc, 
is restored at a'^ (although, unhke other cases seen in the literature [211 ES], the phase 
above a'^ is not exactly the same as that below a^, because of C3 symmetry breaking). 

Analogously to what occurs in a condensation or ferromagnetic phase transition [23] , 
the freezing of the relative occupancy of each state does not imply that individual 
oscillators freeze as well. The frequencies of individual oscillators do decrease with 
increasing a, but only vanish in the limit a — )■ 00, when one of the states is fully 
occupied. 

It is convenient to define an order parameter that is identically zero (in the infinite- 
size limit) for a > a'^. Since the angular frequency u does not fulfill the usual 
requirements for an order parameter [21], we propose 



1^1 = — 



N 



(5) 



x=l 

where 5 is the Kronecker delta and 72,. = gj^j^+i is the transition rate at site x (see 
Eq. [1]). Thus {ipl reflects not only the configuration, but the rate at which the latter 
evolves. 

On the complete graph, is the same for all sites x in the same state j. Denoting 
this rate by jj, the order parameter can be rewritten for mean-field (MF) analysis as 



1 2 MF 



.j=o 



- Po7o-Pi7i - P1I1P2I2 - -P272-P070 • (6) 



Note that both the disordered phase (a < a^) as well as the frozen phase (a > a^) have 
stable fixed points, so Pj'jj = P*_^_i'jj+i from Eq. [2l This in turn renders |'?/'| = in 
Eq. [6] for both cases [a similar line of reasoning can be applied directly to Eq. ([5])]. 

Figure m shows the variation of (= {{\ip\)t)s) with a in mean-field theory, and 
on the complete graph, confirming that functions as an order parameter to detect 
both phase transitions. The MF critical behavior is (lipl) ~ (a — a^Y^"^ for a y/ Uc 
and (l^/'l) ~ {a^ — aY^"^ for a a^. In the pair approximation for a lattice with 
coordination number 2; = 6 we find Oc = 1.97579 and a'^ = 4.71727; for z = 8 the 
corresponding values are Oc = 1.80198 and a'^ = 4.10028. In the pair approximation, no 
transition is observed for z < A. 

On the complete graph, the order parameter decays with system size as {{ipl) ~ 
iV~^/^ at a = ttc and as jsf-^-^'^^^i^) at a = a'^, as shown in Fig. [51 The first result is 
typical of mean-field-like scaling with system size at a continuous phase transition, and 
can be understood as follows. Our numerical results show that with increasing N, 



CONTENTS 



8 



0.25 



0.2 



0.15 



0.1 



0.05 





1 1 


N=2000 






N=5000 






N=10000 






N=20000 


- 




MFA - 




ft 

% 

/••I 

i 1 


\ 

% 

1 1 ' 



1 1.5 2 2.5 3 3.5 

a 

Figure 4. Order parameter as a function of coupling a in mean- field theory and 
on the complete graph, for sizes as indicated. 



as expected, the stationary probability distribution in the Pi — P2 plane becomes 
(for a = ttc) increasingly concentrated around the symmetric point (1/3,1/3). The 
maximum of the distribution, Pmax, scales as N'^^"^ and the number of points at which 
the probability is relatively large (i.e., p{Ni, N2) > Pmaxf^) grows ~ N'^^'^. In terms 
of the densities Pi = Ni/N, the probability distribution is concentrated on a region of 
radius p ~ N~^^^. Note that \-ip\'^ is a quadratic form of the variables = Pi — 1/3 and 
that ip = for ^1 = ^2 = 0. Thus for the purposes of scaling analysis we may write the 
stationary average of {ipl as 

m)- ^%-J -N-^/' (7) 

with ^ = ^/ C,i + ^2- the first transition (a = ac) the same argument applies to the 
order parameter r, and we indeed find numerically r ~ A^""'^^. 

At the second transition, the probability distribution is peaked near the three stable 
fixed points of the mean-field analysis, and so has a very different form than that found 
at Qc- It is therefore not surprising that {\ip\) decays with a different exponent. The value 
of 0.4203(3) (which does not correspond to any simple ratio of integers) nevertheless 
remains something of a puzzle. We suspect that understanding this result will require an 
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Figure 5. (|-!/'|} at the critical points versus system size N on the complete graph. The 
order parameter decays as a power law of the system size at both phase transitions. 



asymptotic (large- iV) analysis of the Fokker-Planck equation for the probability density 
in the Pi — P2 plane, a task we defer to future work. 

3.2. Hypercubic lattices 

Motivated by the appearance of a novel phase transition in mean-field theory and on the 
complete graph, we search for such a transition on finite-dimensional structures, i.e., the 
simple cubic lattice and its four- dimensional (hypercubic) analog, in systems of N = L'^ 
sites. On these lattices, we do observe a marked tendency for u to fall rapidly with a for 
values well above the first transition at (etc — 2.345, ~ 1.900 and ~ 1.750 for d = 3, 
4 and 5, respectively [12] )• We also observe a broad maximum in Xr at some point in 
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Figure 6. (Color online) Four dimensional hypercubic lattice: Inset mean period T 
versus a for (upper to lower) L — 28, 32, 40, and 50. Main graph: scale period T* = 
L'^T; system size increasing upwards. 



this range. To infer the existence of a phase transition, rather than a rapid but smooth 
variation of u and other properties with a, we require evidence of emergent singularities 
as the hnear system size L is increased. 

The global dynamics may be characterized in terms of a process that begins with all 
sites in state j, and ends when the fraction of sites remaining this state falls to < 1/3. Let 
the "escape time" T denote the mean duration this process, averaged over realizations 
of the stochastic dynamics. If there is an infinite-period phase transition at some critical 
value a'^, then T should grow ~ at this point, where z denotes the dynamic critical 
exponent. Although simulations reveal an exponential growth of T with a at fixed L, 
they show that for large, fixed a, the escape time decreases with increasing system size! 
The mean period (see Fig. |6]) exhibits two regimes for a > Qc- for smaller a, relatively 
slow, exponential growth with T e'^ essentially independent of system size, and, for 
larger a, rapid exponential growth with T ~ L~'^e'^"', with k ^ 1. We identify the latter 
regime as that of single-droplet nucleation [25]. In the former regime, the density of 
"advanced" sites (i.e., in state j + 1 in a configuration with the majority in state j) is 
large enough that essentially no barrier to the cyclic dynamics exists. The (smooth) 
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100 




Figure 7. Square lattice, L = 10: evolution of A'^i (main graph) and (inset) in a 
typical realization with a = 8. 



crossover between these regimes occurs at a larger a value for larger system sizes. 

The results for the mean transition rate T on the four-dimensional hypercubic 
lattice, T ~ L"'^ exp[/ta], clearly suggest a nucleation mechanism. The exponential 
factor represents a barrier to nucleation, that is, the mean time to formation of a critical 
cluster. (By a "critical" cluster we mean one that is equally likely to grow as to diminish 
in size; larger clusters tend to grow while smaller ones tend to shrink.) Consider a region 
in which all spins are in the same state. The rate at which a given spin flips to the next 
state (e.g., 1 — > 2) is e~" leading to the exponential growth in T. The prefactor 
corresponds to essentially independent contributions due to each site in the system. 
(Given the extremely low density of clusters in this regime, we can treat the nucleation 
and initial growth of a cluster as occurring in isolation, without interference from other 
clusters.) 

The constant k in the exponential reflects the probability of formation of a critical 
cluster, which in turn depends on the critical cluster size, i*. Intuitively, the nucleation 
barrier exists because an isolated advanced site has a higher rate to advance once again, 
and then rapidly return to the majority state, than it does of inducing a neighbor to 
advance. Given the difficulty of visualizing four-dimensional space, we begin by studying 
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nucleation on the square lattice. Of course, the model does not exhibit an ordered phase 
in two dimensions, but this issue can be bypassed by using small systems (a lattice of 
10x10 sites, in the present case). Then the system remains near one of the three 
attractors (almost all sites in the same state) except for the rapid transitional periods 
following a nucleation event. (In larger systems global synchronization is quickly lost 
and the system breaks into a number of domains.) On such a small lattice, one may 
easily identify transitional regimes between, for example, states with majority-1 and 
majority-2. We define the transition (in the macroscopic sense) as the moment when 
20% of sites are in state 1, with nearly all the rest (70% or more) in state 2. To estimate 
T we determine the time required for a large number (4000, in practice) of transitions to 
occur, in event-driven Monte Carlo simulations. A typical evolution is shown in Fig. [71 
We find that the mean transition time grows exponentially with a; for a > 7 we find 
n = 1.68(1). This suggests that the critical cluster size is approximately 3-4. 

Several other pieces of evidence support this estimate. First, in the pre-transition 
regime, the number N2 of sites in state 2 is seen to fluctuate between zero and three 
or so. In the example shown in Fig. [71 the second time N2 grows to four, a transition 
occurs. Typical formats of a growing cluster are shown in Fig. [HI the cluster is seen to 
be compact, as expected for a growth probability that is highest for the sites outside 
the cluster having the largest number of neighbors within it. 




Figure 8. (Color online) Typical evolution of a cluster of 2s (red) in a background of 
Is (blue). Square lattice, i = 10, a = 6, times t = 39 (left) and 58 (right). 

A semiquantitative analysis of the cluster size distribution in the pretransition 
regime can be developed as follows. Consider first the density p2 of isolated sites in 
state 2, in a background of Is. The rate for a given 1 to change to 2 is e~", while an 
isolated 2 (with all neighbors in state 1) flips to state 3 at unit rate. There are also 
transitions of the form 21 — )■ 22. The total rate for an isolated 2 to become a dimer is 
Thus we have a rate equation of the form 

P2 = e-^pi - P2 - 4e-'^/ViP2 ^ e-" - (8) 

leading to p2 — e^" in the quasi- stationary (QS) regime. By the QS regime we mean the 
long-lived metastable state in which the great majority of sites are in state 1. Similarly, 
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let P22 denote the probability that a given site is the leftmost (or lowermost) site of an 
isolated dimer (the restriction is to avoid double counting). Then we have 

P22 ~ Ae--'^p2 - 2e-^'^p22 - 0{e''''^p22) (9) 

so that P22 — 2e~^"/'^ in the QS state. Consider next the density p\_ of L-shaped trimers. 
The gain term in the equation for p\_ is 4e~'*/^p22 while the dominant loss term (due to 
formation of a four-site, square cluster) is simply —p\_ (the transition rate is unity). Thus 
Pl ~ 8e~'''^/^. The exponential factor, e"''"/^, is quite close to the factor e~^-^^(^)" found 
in simulations, suggesting that this is essentially the critical cluster. (Note that this is 
the smallest cluster with a larger probability of growing than of shrinking.) Applied to 
the simple cubic lattice, this line of reasoning suggests that the critical cluster consists 
of seven sites (in the form of a cube with one corner removed) and yields K3 = 11/3. 
[A crude estimate for the value in four dimensions can be obtained by supposing that 
/t4 ^ {1^3/1^2) X ^^3; this yields k,^ ^ 7.7, reasonably close to the simulation value, 
= 8.43(7).] 

The nucleation scenario outlined above essentially rules out the second transition on 
any finite- dimensional lattice with short-range interactions. There is always a cluster 
size such that growth is more likely than shrinkage, and once such a cluster forms, 
a transition of the global state is likely. From the perspective of equilibrium phase 
transitions, this is not surprising: the ordered state of any (finite) Ising or Potts model 
is always subject to a reversal of magnetization at a finite rate. (This rate, however, 
decays with system size in equilibrium, as a finite cluster is always more likely to shrink 
than to grow, in the absence of an external field.) In the present cyclic dynamics, by 
contrast, there is no impediment to growth (no "surface tension", as it were), once a 
cluster has attained the critical size. 

4. Conclusions 

We have identified a second phase transition in the cyclic model proposed by Wood et 
al. fi2[ \T3\ [m [15] This continuous transition is of a fundamentally different nature 
than the first, as it corresponds to an infinite-period transition at which macroscopic 
properties freeze and a discrete rotational (C3) symmetry is spontaneously broken, in 
the absence of an absorbing state. We find that an appropriate order parameter for 
the second transition involves the local transition rates. While the phase transition is 
observed on the complete graph, we argue that it does not exist on finite-dimensional 
structures with short-range interactions. The question naturally arises whether such 
an infinite-period transition exists in systems with long-range interactions. This indeed 
appears likely, for interactions that decay sufficiently slowly with distance. In the case 
of graphs with interactions between many pairs of sites, independent of distance, we 
expect to observe an infinite-period transition if (in the infinite-size limit) a given site 
interacts with a finite fraction of the others. We plan to investigate these questions, as 
well as the behavior on scale-free networks, in future work. 
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